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Abstract 

Particle-in-Cell (PIC) simulation is the most important numerical tool in plasma physics. How¬ 
ever, its long-term accuracy has not been established. To overcome this difficulty, we developed 
a canonical symplectic PIC method for the Vlasov-Maxwell system by discretizing its canonical 
Poisson bracket. A fast local algorithm to solve the symplectic implicit time advance is discovered 
without root searching or global matrix inversion, enabling applications of the proposed method 
to very large-scale plasma simulations with many, e.g., 10®, degrees of freedom. The long-term 
accuracy and fidelity of the algorithm enables us to numerically confirm Mouhot and Villani’s the¬ 
ory and conjecture on nonlinear Landau damping over several orders of magnitude using the PIC 
method, and to calculate the nonlinear evolution of the reflectivity during the mode conversion 
process from extraordinary waves to Bernstein waves. 

PACS numbers: 52.65.Rr, 52.25.Dg 
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In modern plasma physics, numerically solving the Vlasov-Maxwell (VMl^equations using 


the Particle-In-Cell (PIC) method has become the most important tool [l|, Q] for theoretical 
studies in the last half century. Many innovative algorithms, such as the Boris scheme for 
advancing particles 0,3 and Villasenor-Buneman’s charge-conserving deposition scheme fl, 
have been developed and successfully applied. Recently, new geometric numerical method¬ 
ology has been adopted for PIC simulations. This exciting trend begins with the discovery 
of symplectic algorithms for Hamiltonian equations that govern charged particle dynamics 


6l-ll6l|. The Boris algorithm was discovered to be volume-preserving [17|, [iSj] and high-order 
volume-preserving methods have been found 1^. In addition, the Vlasov-Maxwell sys¬ 
tem 20I424| and the Vlasov-Poisson system 2^ have been discretized from a variational 


symplectic perspective that preserves symplectic structures and exhibits excellent long-term 
accuracy and hdelity. 

In this letter, we develop a new canonical symplectic PIC method for solving the VM 


equations by discretizing its canonical Poisson bracket 2^. The distribution function / is 


hrst discretized in phase space through the Klimontovich representation by a hnite number 
of Lagrangian sampling points (X*, Pj) {i = 1,..., N), where Xj and are the position and 
canonical momentum of the i-th particle, and N is the total number of sampling points. The 
electromagnetic held is discretized point-wise on a given spatial grid, and the Hamiltonian 
functional is expressed as a function of the sampling points and the discretized electro¬ 
magnetic held. This procedure generates a hnite-dimensional Hamiltonian system with a 
canonical symplectic structure. The number of degrees of freedom for the discrete system is 
D = 3N + 3M, where M denotes the total number of the discrete grid-points. 

In general, for a Hamiltonian function whose momentum dependence and position de¬ 
pendence are not separable, it is not possible to make symplectic integration algorithms 


explicit 


0 , 


14i\ . For the discrete Hamiltonian system developed here for the VM equations. 


the dimension of system is usually very large, and root searching algorithms required by 
implicit methods are too time-consuming to be practical. However, we discovered that if the 
symplectic Euler algorithm 9j is applied to the discrete VM Hamiltonian system at hand, 
the implicit time advance can be carried out as inexpensively as an explicit method by just 
inverting a 3 x 3 matrix for every particle separately. The resulting canonical symplectic PIC 
method for the VM system inherits all the good numerical features of canonical symplec¬ 
tic algorithms, such as the long-term bound on energy-momentum error. Being symplectic 


2 






















means that the numerical solution satishes D{2D — 1) constraints as the exact solution does. 
Since D is a large number, the symplectic condition is much stronger than a few constraints 
on global energy and momentum. The symplectic condition is almost as strong as imposing 
local conservation everywhere in phase space. 

Two examples of application are given. In the hrst example, we simulate the dynam¬ 
ics of nonlinear Landau damping. It also serves as a test of the algorithm. The discrete 
VM Hamiltonian system for this study has more than 2.69 x 10® degrees of freedom. The 
damping rate from the numerical results agrees exactly with the theoretical result. Fur¬ 
thermore, long-term simulations reveal that the phase mixing dynamics in velocity space is 


the physical mechanism of the nonlinear Landau damping, as recently proved 


and Villani 


jy Mouhot 


271. l28l | for the Vlasov-Poisson system and conjectured by Villani 


29] for the 


Vlasov-Maxwell system. In the second application, we study the nonlinear mode conversion 
process from extraordinary modes to Bernstein modes (X-B mode conversion) in an inhomo¬ 
geneous hot plasma. Simulations show that nonlinear mode excitations and self-interaction 
of the Bernstein waves signihcantly modify the reflectivity and conversion rate. It is the 
long-term accuracy and hdelity of the canonical symplectic PIC algorithm that enables us 
to numerically conhrm Mouhot and Villani’s theory and conjecture over several orders of 
magnitude using the PIC method, and to calculate the nonlinear evolution of the reflectivity 
during the X-B mode conversion. 


We start 


equations |26|, 


rom the canonical Poission bracket and Hamiltonian for the Vlasov-Maxwell 


{F,G}^j /|-.-|..„dxdp + y 


^?(/,A,Y) = - 


(p - Affdxdp +- + (V X A) 


dx. 


( 1 ) 

( 2 ) 


Here, F, G, and the Hamiltonian H are functionals of the distribution function /, vector 
potential A, and Y = dA/dt. The bracket {h,5f}xp inside the hrst term on the right hand 
side of Eq. ([T]) is the canonical Poisson bracket for functions h and g of canonical phase space 
(x, p). The temporal gauge, i.e., 0 = 0, has been explicitly chosen for this Poisson bracket 
to be valid. The Poisson bracket dehned in Eq. ([T]) can be formally derived from the point of 
view of co-adjoint orbit theory 26(|, and can be used to derive the non-canonical Morrison- 


Marsden-Weinstein bracket in the (/, E, B) and (x, v) coordinates 


2 ^, 


30, 


3^. First, we 
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discretize the distribution function using the Klimontovich representation 

N 

/(x,p,t) = ^5(x-Xi)(5(p-Pi), (3) 

i=l 

where (Xj, Pj) {i = 1, ...,N) are particles’ coordinates in phase space. Under this discretiza¬ 
tion, it can be shown that 

^ = y'nx-x<)np-Pi)^(|^)<ixdp. (4) 

from which we obtain 


6F 6G 6G 6F 


5X, 6Pi 


= j (5(x-X,)5(p-P, 


5F 5G' 


6f' 6f 


5:s.i 5P, 

It then follows that the hrst term on the right-hand side of Eq. ([T]) is 


xp 


dxdp. 


/ 




6G 6F 


2 = 1 


6Xi 5P,: (5X,- (5P, 


( 6 ) 


(7) 


Similar derivation of the discretized bracket can be found in the context of Hamiltonian 
description of vortex fluid 
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To discretize the second term on the right-hand side of Eq. ([T]), we hrst discretize the 
helds A(t) and Y(t) on a Eulerian spatial grid as 


M M 

A(x,t) = ^ Aj(t)T(x-xj), Y(x,t) = ^ Yj(t)T(x-xj), 

,7=1 J=1 


( 8 ) 


where the discrete helds Aj(t) and Yj(t) are the helds evaluated on the grid-point xj. The 
subscript J is the index of the grid-point, and M is the total number of the grid-points. 
Here, T(x — xj) is the step function. 


1, lx - xj| < ^,\y- yj\ < ^,\z- zj\ < ^ , 


T(x - xj) = 


0, elsewhere . 

Under this discretization of A{t) and Y(f), we have 


(9) 


SF " SAj dF 1 ,, , dF 

= E = E —'l'(x - X,)- 


( 10 ) 


SA 5A dAj ^ ' dAj 

where AU is the volume of each cell, which is taken to be a constant in the present study. 
For the second equal sign in Eq. fflll]) . use has been made of the fact that 

1 


5A AU 


T(x-xj). 


( 11 ) 
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The discretization of the second term on the right-hand side of Eq. ([T]) is thus 


r /6F 6G 

J 


SG5F\ ( dF dG 

JaW) ~ [dAj dYj 


dG dF\ 1 


( 12 ) 


Finally, the discrete Poisson bracket for the VM system is 



6F 6G 


5G 5F\ ^ f dF dG dG dF \ 1 

^ Uaj dYj ~ dAj dYj) W 


(13) 


for functions F and G of the particles (Xj,Pj) and the discretized field (Aj, Yj). 

Next, we need to express the Hamiltonian functional given by Eq. ([ 2 ]) in terms of (Xj, Pj) 
and (Aj, Yj). The particles’ total kinetic energy is the sum of each particle’s kinetic energy. 
The vector potential at a particle’s position can be interpolated from Aj{t) as 


M 

A(X„ t) = Y. Aj{t)W{y., - xj), ( 14 ) 

J=1 

where W (X^—xj) is a chosen interpolation function. Note that W (X^—xj) is not necessarily 
the same as the step function T(x — xj) in Eq. ([H]). This is of course allowed as long as the 
consistency condition is satished, i.e., the continuous limit is recovered when the grid-size 
goes to zero. The Hamiltonian then becomes 


N 


M 


H(X„ P„ Aj, Yj) = - j: ^ p2 - 2P, ■ ^ AjPE(X, - xj) 


i=l 


M 


M 


^ Aj ■ AlW{X, - xj)fE(X, - Xi) I + - ^ [y2 + (Vj X A)( 
J,L=1 J ^ J=1 


AV, (15) 


where (Vj x A)j is the discrete curl operator acting on the discrete vector potential eval¬ 
uated at the J-th grid-point. Finally, the discrete Hamiltonian flTHll and discrete Poisson 
structure 0131) form a canonical symplectic discretization of the original continuous Vlasov- 
Maxwell system. The ordinary differential equations for the canonical system are 
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M 


X, = 

{Xi, H\ = p. - E AjfE(X, - xj), 

(16) 


j=i 


A, 7 = 


(17) 


M 


P* = 

{P- = E (P* ■ Aj) V1E(X, - xj) - 



J=l 


M 

(Aj-Ai)iy(X,-xj)Viy(X,-x^), (18) 

J,L=1 

^ 1 

Y^ = {Y./.J^}, = gPWX.-x.)^- 

N M ^ 

Y, E AiH'(X, - xj)iy(X. - Xi)— - (vj X Vj X a) (19) 

i=lL=l IXV ^ 'J 


This equation system consists of 6(M + A) equations describing the dynamics of N particles 
and fields on M discrete grid-points. The last term in Eq. is dehned to be 


(Vj X Vd X A)^ 


1 d 
2M 


- M 

E (Vd X A) 


.L=l 


2 

L 


( 20 ) 


The notation ofVjxV^xA indicates that the right hand side of Eq. fl2UI) can be viewed 
as the discretized V x V x A for a well-chosen discrete curl operator V^. To wit, we note 
that the term V x A in the Hamiltonian is discretized using the step function \h(x — xj), 


M 

Vx A = 5](VrfX A)^T(x-xj). 


j=i 


( 21 ) 


As an example, we dehne the discrete curl operator (V^ x A)j = (V x A)^^. ^ to be 


(V. X A), ^ 


A .. , , 

t,j,k — 

2^^ • ■ u 1 

2,7,fe 2,7,fe —1 

Ay 

Az 

1 

1 

'E? 

1 

Az 

Ax 

^ t- 

i,7,fe 2—l,7,fe 

1 

'E? 

1 

Ax 

Ay 


( 22 ) 


which can be written as a linear operator on the space of 3M-vectors as 



^ (V, X A)i ^ 


^ A, \ 

Vd X A = 


= T 



(Vd X A.)j^ j 


V Am j 


( 23 ) 
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Here, F is a 3M x 3M sparse matrix specifying the discrete curl operator. The partial 
derivative with respect to Aj can be expressed as 

/a, \ 


1 d 


2dAj 


M 


yvjx A)i 

.L=l 


r^r 


\ Am / 


“ (Vl X Vd X a)^, 


(24) 


where F^ is the transposition of F. Obviously, the notation in Eq. fl2U]) or fl2T)) is meaningful 
for any linear discrete curl operator V^. 

It is clear from Eqs. fll6l) - ffT^ that the particles and helds interact through the interpola¬ 


tion function hF(Xj — xj). The function hF(Xi — xl)/AV distri 


the grid-points as if they are “charged clouds” with hnite-size [l|. 

Once the canonical symplectic structure is given, canonical symplectic algorithms can be 
readily constructed using well-developed methods |6l-ll3j|. For a reason soon to be clear, we 
adopt the semi-explicit symplectic Euler method for time advance. The symplectic Euler 
method for a generic canonical Hamiltonian system is 

dH 


Dutes particles’ charge over 


= p- 


At 


dq 

dH 




= g’' + At— 

op 


(26) 

(26) 


where At is the time-step, and the superscript n in p” and g” denotes that they are the value 
at the n-th time step. It is implicit for p, but explicit for q. Making use of this algorithm, 
the iteration rules for Eqs. (Il6|) - (Il9l) are 


vn+l _ vn M 

— —^=py - E V3»'(x” - xj), 


At 

A«+l — \ri 

_ -vx^+l 

At ’ 

pn+l _ pn M 


J=1 


(27) 

(28) 


At 


E (PX' • a;) vh'(x” - X,,) 


J=i 


M 


- E (AV a;) W'(X” - X,,)VH/(X" - xi), 

J,L=1 

N 1 

•'= j;pyw/(X”-x,,) 


(29) 


At 


2 = 1 


AH 


N M 


-EE A2M-'(X“ - xj)»'(X” - Xi)^ - (vj X V, X A 


i=l L=1 


(30) 
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These difference equations furnish a canonical symplectic PIC method for the Vlasov- 
Maxwell equations. 


As discussed above, symplectic algorithms for a Hamiltonian system with non-separable 
momentum and position dependence are implicit in general. This is indeed the case for 
the difference equations (l27D - fl5(l . because the right-hand sides of Eqs. depend on 

values of the (n -f l)-th time-step. However, they are semi-explicit, because Eqs. (ET]), (l28jl . 
and are explicit for and respectively. Another good property of the 

system is that the only implicit equation 0291) is linear in terms of and it is only implicit 

for each particle, i.e., Eq. does not couple and P^"'"^ when i ^ k. Therefore, the 

system can be solved without root searching iterations as follows. We first solve the linear 
equation (12^ for P^^^ for every index i separately, which amounts to inverting a 3 x 3 
matrix for every i. Then and Yj"*"^ are advanced explicitly according to Eqs. 0271) and 
(j3n|) . and the last step is to advance also explicitly, according to Eq. (j28i) . 


The preservation of the symplectic structure exerts D{2D—1) constraints on the numerical 
solution. Because D is a large number, preservation of symplectic structure is a very strong 
constraint and signihcantly reduces the errors of numerical solutions. We can also appreciate 
this advantage from the viewpoint of symplectic capacity, which is defined on any open set of 


the phase space. Symplectic maps preserves symplectic capacity 


34l |. and in principle there 


are inhnite constraints that symplectic algorithms can satisfy as the continuous systems do. 


We now apply this canonical symplectic PIC scheme to simulate the nonlinear Landau 
damping process. This study also serves as a test of the algorithm. Previously, similar 
study and test have been performed for other algorithms, e.g., the Eulerian algorithms for 
the Vlasov-Poisson system 22|, [3^. The ions are treated as a uniform positively charged 
background, and the dynamics of electrons are simulated. The electron density is Ug = 
1.149 X 10^®/m^, and the thermal velocity of electrons is vt = 0.1c, where c is the light 
velocity in vacuum. The three-dimensional computational region is divided into 896 x 1 x 1 
cubic cells. The size of grid is chosen to be Ax = 2.4355 x 10“^m, the time step At = Ax/2c. 


8 






The interpolation function is chosen to be 8th order, i.e., 

lT(x) = Wi{x/Ax)Wi{y//\x)Wi{z/Ax) , 

0 , 


( 31 ) 


n\{q) = 


15 8 _ 1^„7 _ ^„5 ^„4 _ ^ 1 

^ 128'^ “T 19sy ' (\A^ y “T 


1024 


128 


32 ■ 


64 


q>2 , 
l<q<2, 


15 ^8 


0» _ T^o7 , _Lo 6 _ ^ 0 , ilO 4 _ iUO 2 , iil 0 < 0 < 1 

_y 128^ ^ 16^ 32^ ^ 256^ 128^ ^ 512’ ^ y ^ J- , 

--|- T^o'^ + —O® + —Q^ + —n'^ — —Q^ + — —1 < a < 0 

1024y ^ 128^ ^ 16^ ^ 32'J ^ 256^ 128^ ^512’ ^ y ^ ) 

-I- _j_ _4^^6 _j_ ^^5 _j_ _|_ ^ _|_ 1 —9 n ^ —1 

1024 " ^ 128 128 32 "^ 64"^"^'^’ ^ q — Ij 

0, q<-2 . 

It can be proved that the kernel function W is 3rd order continuous in the whole space. 
According to our performance benchmarks, this 8th order kernel is about 30% more compu¬ 
tationally costly than a 2nd order kernel, which is acceptable since a higher order continuous 
kernel gives more numerical fidelity. Initially, 10^ sampling points of electrons are distributed 
in each cell. The total number of particles is = 8.96 x 10^, and the number of degrees of 
freedom is D = 2.69 x 10®. The initial electrical field perturbation is Ei = Ei cos {kx) e^, 
where the wave number is k = 27r/224Aa;, and the amplitude of the perturbation electric 
held is El = 9.103 x lO^ld/m. The simulations are performed for 80000 time steps, dur¬ 
ing which a complete picture of the nonlinear Landau damping is revealed. As expected 
for symplectic algorithms, the numerical error of energy does not increase with time and is 
bounded within 1% for all time. The theoretical damping rate calculated from the dispersion 
relation is Ui = —1.3926 x 10^/s, and the theoretical real frequency is Ur = 9.116 x 10®/s. In 
Fig. [U the slope of the green line is the theoretical damping rate, and the blue curve is the 
evolution of the electrical held observed in the simulation. It is evident that the simulation 
and theory agree perfectly. After t = 30/0;^, the energy of the wave drops below the level of 
numerical noise, and the damping process stops. The evolution of the electron distribution 
function is plotted in Fig. [2|, which clearly demonstrates the mechanism of phase mixing 
in velocity space. We observe in Fig. [2] that the wave-number in velocity space increases 
with time, which results in a decrease in density perturbation and thus attenuation of the 
electrical held. More importantly, this mechanism of phase mixing is the dominant physics 
or the entire nonlinear evolution of the Landau damping, as proved by Mouhot and Villani 
27 . ^ recently for the electrostatic Vlasov-Poisson system. In addition, our simulation is 
electromagnetic and it shows that this physical picture of nonlinear Landau damping is also 


7 ^6 21 ^5 


175 ^4 


105 ^2 


337 


(32) 


valid for the electromagnetic Vlasov-Maxwell system, as Villani conjectured 


29|. It is the 
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FIG. 1. Perturbed electrical field as a function of time. The slope of the green line is the theoretical 
damping rate. 



FIG. 2. Electron distribution function in velocity space as a function of time. Different colors 
denote the amplitude of the perturbation. The mechanism of phase mixing in velocity space is 
clearly demonstrated. The wave-number in velocity space increases with time, which results in a 
decrease in density perturbation and thus attenuation of the electrical field. 

long-term accuracy and fidelity of the canonical symplectic PIC algorithm that enables the 
confirmation of Mouhot and Villani’s theory and conjecture over several orders of magnitude. 

Even though Mouhot and Villani rigorously proved that when the initial perturbation 
amplitude is small enough, the electrical perturbation will decay to zero, plasma physicists 
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FIG. 3. In nonlinear Landau damping, the perturbation will bounce back after initial phase of 
damping if the initial perturbation is large enough. The amplitude of the initial electric field in 
this case is 0.494MV/m. 


have known this fact since 1960s 


36l. l37( |. It has also been known that when the initial 


electrical perturbation is large enough, the perturbation will bounce back after initial phase 
of damping [36|,l37t]. One such case simulated is plotted in Fig.[3l For this case, the amplitude 
of the initial electrical held is 0.494MV/m and the wave number is k = 27^/2'12/S.x. Figure 0] 
shows the bounce-time as a function of the initial amplitude of the electrical held obtained 
in simulations. The physics demonstrated in our PIC simulations agrees with that obtained 


35[ |. except that our simulations are carried out for the full Vlasov- 


from Eulerian solvers 
Maxwell system. 

In the second application, the ID nonlinear mode conversion of extraordinary waves to 
electron Bernstein waves (X-B mode conversion) in a inhomogeneous hot plasma is simu¬ 
lated for a long time. The plasma density prohle in the simulation is 


ne(x) = no 



^ ^ 320) , 

K + 320) <f^< 1300 , 


(33) 


where no = 2.3 x 10^® m“^, n^ = 380, and Ax = 2.773 x 10“^m is the grid size. The 
thickness of the plasma boundary is rirAx. The electron temperature is Tg = 57.6eV, and 
constant external magnetic held is B = BqBz with Bq = 0.6T. The simulation domain 
is a 1584 x 1 x 1 cubic mesh. At both boundaries in the x-direction, the Mur’s absorbing 
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FIG. 4. Relationship between the initial amplitude Ei and the bounce time tb in nonlinear Landau 
damping. 


condition are nsed, and periodic bonndary conditions are adopted in the y- and 2 ;-directions. 
The time step is chosen to be At = Ax/2c. At the left bonndary, a sonrce is placed to excite 
an electromagnetic perturbation at a; = 0.0145/At with Ei = Eiey and Ei = OOOkV. As 
illustrated in Fig. El the extraordinary wave excited at the left boundary hrst propagates to 
the region of cutoff-resonance-cutoff near xjAx = 500. The wave is then partially reflected 


back, and partially converted to electron Bernstein waves 


24l |. The reflectivity evolution 


is plotted in Fig. O Nonlinear excitations and self-interactions of the Bernstein modes 
dominate the long time dynamics of the mode conversion process. As a consequence, the 
reflectivity of the incident wave evolves nonlinearly as shown in Fig.O For this set of chosen 
parameters, the incident wave is completely reflected at the later time of the process. 

In conclusion, we have developed a canonical symplectic particle-in-cell simulation 
method for the Vlasov-Maxwell system by discretizing its canonical Poisson bracket. In 
phase space, the distribution function is discretized by the Klimontovich representation us¬ 
ing Lagrangian markers, and the electromagnetic held is discretized point-wise on a spatial 
grid. The resulting canonical Hamiltonian system with a large number of degrees of free¬ 
dom is integrated by the symplectic Euler method, whose difference equations can be solved 
inexpensively by inverting a 3 x 3 matrix locally for every particle. Implicit root searching 
and global matrix inversion are avoided entirely. This technique makes large-scale applica¬ 
tions of the developed canonical symplectic PIC method possible. To suppress numerical 
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FIG. 5. The space-time dependence of Ey(t,x) during the nonlinear X-B mode conversion. 



FIG. 6. The evolution of reflectivity during the nonlinear X-B mode conversion. 


noise caused by the coarse sampling, smoothing functions for sampling points can also be 
conveniently implemented in the canonical symplectic PIC algorithm. By incorporating the 
smoothing functions into the Hamiltonian functional before the discretization, we are able to 
rein in all the benehts of smoothing functions without destroying the canonical symplectic 
structure. Progress in this and other directions will be reported in future publications. 
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